Synchronization of oscillators with long range interaction: phase transition and 

anomalous finite size eff'ects 



Mate Marodi,^ Francesco d'OvidioJ^ and Tamas Vicsek|] 
- - - *XDept. of Biological Physics, Eotvos University, 

(N ■ Budapest, Pdzmdny P. Stny. lA, 1117 Hungary 

, ^ Chaos Group and Quantum Protein Center, Dept. of Physics, 

>^ ' Building 309, Danish Technical University, DK-2800 Lyngby, Denmark 

^ (Dated: February 1, 2008) 

] 

^ Synchronization in a lattice of a finite population of phase oscillators with algebraically decaying, 

non-normalized coupling is studied by numerical simulations. A critical level of decay is found, below 
which full locking takes place if the population contains a sufficiently large number of elements. 
OA ' For large number of oscillators and small coupling constant, numerical simulations and analytical 

OA , arguments indicate that a phase transition separating synchronization from incoherence appears at 

a decay exponent value equal to the number of dimensions of the lattice. In contrast with earlier 
results on similar systems with normalized coupling, we have indication that for the decay exponent 
^ , less than the dimensions of the lattice and for large populations, synchronization is possible even if 

. the coupling is arbitarily weak. This finding suggests that in organisms interacting through slowly 

decaying signals like light or sound, collective oscillations can always be established if the population 
is sufficiently large. 



(N 



< 



a 

o 
o 



> 



O 



X 



PACS numbers: 05.45.Xt 



I. INTRODUCTION 



Synchronization is one of the most fascinating nonlinear phenomena appearing in a wide range of fields. Examples 
cover physical systems (network of Josephson junctions [Q), oscillating chemical reactions molecular turnover in 
allosteric enzymes Q , a variety of biological observations (synchronous flashing of fireflies Q , menstrual synchrony , 
metabolic activity in yeast cells ^, and other various physiological processes Q), and social phenomena (rhythmic 
applause |^). 

Models for rhythmic synchronization use populations of nonlincarly coupled oscillators [|lO|, |ri|, [T^ . These models 
' might be grouped according to the nature of coupling. For instance, for pulsatile interaction integrate-and-fire 
fT^ , oscillators provide a useful tool of description ||l3whilc in the case of continuous interaction a system of coupled limit 
I ■ cycle oscillators has proved to be a good modeTp]. 
^ An important aspect of synchronizing systems is the spatial dependence of the coupling between two oscillators. 
•■ Two limit cases have been deeply investigated: i) mean field (with each pair of oscillators interacting with a given 
coupling strength, independent of their position, ii) oscillators on a lattice with interactions only between nearest 
neighbours. However, realistic systems are quite different. In many of the biologically relevant systems the signals 
carrying information about the phase of an oscillator decay as a function of the distance slowly. For example, the 
^ , intensity of the sound or light signals decays as a power law, in particular, in the three dimensional case as the inverse 
■ of square of the distance from the source. In Ref. results have been presented concerning the existence of a critical 
^ [ exponent below which full synchronization can be achieved. These results considered a coupling term normalized by 
O ■ the sum of the spatial coefficients and the limit of the population size going to infinity. 

However, in a real system typically the interaction between two oscillators depends only on the spatial distance 
and the phase differences and thus there is no normalization coefficient. This change, trivial at first sight, introduces 
a peculiar dependence of the dynamics on the number of oscillators, and thus places the study of finite populations 
. . into relevant context. 

In this paper synchronization in populations of phase oscillators with a non-normalized coupling term is studied. 
Our investigations are especially focused on the dependence of full locking and the mean field on the number of 
oscillators for different value of the exponent a describing the spatial decay of the interactions. The main result of 
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our work is the identification (through numerical simulations and a simple analytical argument) of critical values 
of a below which full synchronization, or in more general, collective oscillations can take place if the population 
is sufficiently large. We remark that very few works in our knowledge has approached the problem of finite size 
populations. 

The paper is organized as follows. In S ec. II we summarize previous results on the versions of the model of our 
studies, the Kuramoto model. In Sec. Ill we explore the boundary of full locking in the plane of the decay exponent 
a and the number N of oscillators, at different coupling values. We find that if a is below a critical value ac{K) full 
synchronization can be achieved if the system is sufhciently large, i.e., above a critical size N. We argue however that 
this information is partial for real systems, since bulk oscillations can take place also if the system is not completely 
locked. We thus study in Sec. IV the behaviour of the mean field. It appears then that a phenomenon close to a 
thermodynamic phase transition takes place looking at the mean field changing a, with a discontinuity developing for 
N going to infinity and increased fluctuations around the transition point. For K going to zero, analytical arguments 
and numerical simulations indicate that the transition point occurs at a equal to the number of dimensions of the 
lattice. To show the robustness and generality of the result, we finally study a system of phase oscillators where the 
diversity is given by thermal noise. In such system full locking is not strictly defined, but on the contrary a mean field 
approach can be carried out in the same way as for oscillators with fixed parameter distribution and the same kind of 
phase transition appears. The system with thermal noise is very similar to the Heisenberg model with ferromagnetic 
interaction. We shall suggest that a thermodynamic formalism may be applied to study analytically the phenomenon 
discussed numerically in this paper. 

We remark that the present model has been proposed for real populations of oscillators (among others, specific 
examples are the above mentioned works on Josephson junctions ||l| and synchronized applause [^). In these cases, 
and when spatial dependences are included, our results have direct implications. They suggest that whenever the 
coupling signals of these systems decay slower than the number of spatial dimensions (like for example light or sound) , 
collective oscillations always arise if the population is sufficiently large, however weak the coupling is. 



II. THE KURAMOTO MODEL: PREVIOUS RESULTS 

The original form of the Kuramoto model was introduced to describe oscillating chemical reactions ^ . Later 
it proved to be useful in modeling a wide variety of processes (see e.g. Rcf. |l^). The basic concept is a population 
of N interacting oscillators coupled via nonlinear phase-difference minimizing interaction. In general the equation of 
motion of the phase ipi of the i-th oscillator is: 



(1) 



where iOi is the natural frequency of the i-th oscillator, r(0) is the two-oscillator interaction and the sum goes over 
a suitably defined subset of the population depending on the actual model. The system can be further simplified 
choosing, without loss of generality, a frame of reference rotating at the average natural frequency Ui — loq and 
rescaling the time and the coupling constant by the inverse of the standard deviation of the natural frequencies. The 
distribution of natural frequencies becomes then centered in zero and its standard deviation normalized. 
Macroscopic states of this system can be characterized by a real order parameter: 



z(t) 



N 



(2) 



Another approach to synchronization is the investigation of frequency- locked states, i.e. when all 4'i{t) approach 
asymptotically the same value. 

Two basic setups of the Kuramoto model have been studied extensively so far. First, let us consider the case when 



r,,(0)--sin((/.), (3) 

which defines a mean-field interaction^ where K > Q \s the coupling strength. Now the summation in (|l|) goes over 
the whole population. Let us assume that the distribution of the natural frequencies oji is symmetric about 0, and is 
convex in that point. Under these conditions the existence of a critical coupling was proved, above which partial 
synchronization is possible ||l^, |2^. With the above conditions the frequency of the locked subset is equal to the 
mean of the natural frequency distribution. 
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The other well explored case is when the interaction is inherently local, i.e. nearest neighbours are coupled. Then 
the model takes the form: 



sin((/)j - (4) 

where the symbol denotes summing over all the nearest neighbours of the z-th oscillator. Strogatz and MiroUo 

proved that the probability of a phase-locked solution tends to zero as the number of oscillators in the system goes to 



infinity |21 . Also, they studied the more general case when the mean interaction exerted on one oscillator is uniformly 
bounded and the distribution of natural frequencies "sufficiently broad-banded" (for details see Ref . |^ ) . For these 
conditions they showed that in the N oo limit the probability of a phase-locked state is zero. 
An interesting situation arises when the interaction decays in space as a power function: 

=Wi-| ^ — sin (5) 

where is the distance between the i-th and j-th oscillator and 77 is a normalization coefficient. It is clear that Eq. 

is equivalent to the mean-field case when a ~ 0. Also, when a — > 00, the model becomes the same as in the 
nearest neighbours case. So, changing the exponent a, one can have a continuous transition between the two extremes. 
Knowing the result of full synchronization in the global coupling case for high and only local synchronization in the 
nearest neighbours coupling for any K, one may expect a critical value of a below which full locking may be achieved. 
This was shown by Rogers and Wille, who studied Eq. (^) setting 77 = Z^j^j"^-*^^ in the hmit of N going to 
infinity Jl^ . They demonstrated that critical values of a exist depending on the coupling strength K. A similar 
result holds for sine-circle maps as well |p2| . It was also recently shown numerically that an adequate normalization 77 
removes the dependency of the exponent a and system size N on the fraction of oscillators synchronized in frequency 
in the case a < d |p3| . 

However, in most real systems the coupling between two oscillators does not depend on global information on the 
population, but only on the phases and the distance of the two oscillators. It is thus realistic for these cases to set the 
normalization coefficient to 1. As an important consequence, the behaviour of the oscillators then strongly depends 
on the size of the population. 

Throughout this paper we study models where rj = 1. Before describing any numerical results, we remark that an 
estimation of a critical point can be obtained looking at Eq. and considering the coupling term. It is then natural 
to expect that the point a — d, where d is the number of dimensions of the lattice, has to play a special role. In fact 
if a > d the coupling term is bounded, for every N: 



^ 1 



00 ^ 



However, for a < d, in the limit of large N the value of the coupling term may entirely dominate the difference in 
natural frequencies LOi . More precisely, there are regions of the phase space where this term may diverge for N going 
to infinity, in particular, in those regions where the difference between phases is large. Hence, we expect that for a < d 
synchronization is enhanced enlarging the size of the system, independently of the (positive) value of the coupling 
constant K. This property and the effects concerning the population size result from setting the coupling without 
normalization, and follow directly from the divergence of the upper bound in relation Eq. (^. 



III. FINITE POPULATION SIZE AND FULL SYNCHRONIZATION BOUNDARIES 

In this Sec. we explore the effect of changing the number of oscillators in a one-dimensional system, measuring 
the boundary of the full locking region in the (a, N) plane, for different K values. From the brief discussion in the 
previous Sec. we may expect that for sufficiently low values of K synchronization can be achieved if N is sufficiently 
large, only if a < 1. On the contrary, this should not be possible for a > 1. 

In order to identify the exact value of a at which the transition occurs it is helpful to consider the fraction p of the 
oscillators that are locked at the main frequency uq. We thus count how many oscillators asymptotically satisfy the 
condition: 
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(7) 



and we divide this number by N. 

The plane {N, a) has been scanned in either vertical or horizontal direction, numerically integrating population 
of oscillators at those parameter values. The transition points were identified as the boundary of the region of full 
synchronization, i.e., measuring the transient p and finding the parameter values where p is no longer equal to f . 
Refinements up to Aa = O.Of and AiV = 1 were used. In all numerical computations in this paper the natural 
frequencies tOi were chosen from a Gaussian distribution with expectation value oJ,; — and unit variance. Since the 
value of Nc can in principle depend on the position of the natural frequencies on the lattice, each point has been 
averaged on several configurations, especially for low values of N. We remark that, setting without loss of generality 
the main frequency to zero, the locked state appears as a steady state. To locate it, a direct integration method 
(Euler) has been used. This method is simple to implement and reliable, since Eqs. (||) are clearly non-stiff, but 
requires a careful setting of the parameters for computing p: in fact, the line Nc corresponds to bifurcation points, 
where the transient time becomes more than exponentially long. Other methods, like continuation, could give a better 
estimation of the line but may be difficult to apply for large systems, having to deal with matrices of the order of 

For high coupling, full synchronization can appear at every a values. However, if K is low we expect synchronization 
only for a values below the dimensions of the lattice. Results are plotted in Fig. The region a > 1 is outside the 
boundary of full synchronization, confirming the discussion of Sec. I. The situation is more complex for a < 1. Here 
the effect of reducing the coupling can be indeed compensated by increasing the number of oscillators, up to a critical 
value of a. Such critical value however seems to be in general different from I (the number of dimensions of the 
lattice), and depending on the coupling strength K. This is due to the fact that for small K values, synchronization 
is still enhanced by a larger number of oscillators, although the system is not completely locked. The next Sec. will 
investigate this situation, looking at the behaviour of the mean field. 
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FIG. 1: Boundary of the full locking region in the (A'^, a) plane for different values of K. The line a = 1 is also shown. 
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IV. COLLECTIVE OSCILLATIONS 

The results of the previous Sec. may be relevant for real systems. They suggest that in the case of oscillators 
coupled with a slowly decaying signal, full synchronization can be achieved as soon as the population is sufficiently 
large. However, the information that they give, describes only partially what may appear macroscopically. In this 
Sec. we propose to approach the problem looking at the behaviour of the mean field when a is changed. This criterion 
is less sharp than full synchronization: a transition cannot be seen for finite N but, in analogy with thermodynamic 
phase transitions, looking at discontinuities for N going to infinity. Nevertheless, as we shall see, the order parameter 
approach provides a robust and meaningful framework for describing the relations between synchronization and the 
decay of the coupling. Moreover, it allows to state a result in a simple way, connecting the critical value of the 
decaying of interaction with the number of dimensions of the lattice. 

Let us consider numerical simulations conducted on one and two dimensional lattices. We integrated the system 
of coupled differential equations using the Euler method, timesteps varied from 0.01 to 0.0005 relative time-units 
(measured in comparison with the frequencies). The value of the order parameter is computed calculating at each 
time z{t) = |l/iV^^-^j^ e*"^J*^*^|, discarding the first few thousand steps of integrations as transient. Averaging this 
value on a long time (typically, 20000 time steps) gives the measure of synchronization we use, Z ^ {z{t)). 




FIG. 2: Time average of the order parameter as a function of a in d = 1. Simulation results for = 50, 100, 200, 500, 1000 
(symbols are squares, pluses, diamonds, circles, and down triangles respectively). The coupling is if = 1.0. 

One-dimensional results are plotted on Fig. ^, where one can observe the transition from highly synchronized states 
at low a values to unsynchronized states. The behaviour of the system can be divided into three regimes. For a ^ 1 
the system is fully locked, while for a ^ 2.5 the order parameter Z approaches a steady value. In between one finds the 
region of transition. These three regions represent three different microscopic behaviors corresponding to the different 
macroscopic states (see Fig. ||). When a is low (e.g., a = 0.5) the system rapidly reaches the state of complete 
synchronization, when for all oscillators their frequencies become 0j — loq = 0. The next two subfigures represent 
the transitional region {a ~ 1.5,1.8). Here one can observe clusters arising first with increasing phase differences, 
then with synchronization holding only for finite times. This phenomenon leads to large fluctuations in the order 
parameter. Finally, in the unsynchronized region {a = 2.5) there are only local synchronized groups, and clearly only 
very close oscillators influence each other. 
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The transition can be observed also in two dimensions (see Fig. ^). Simulations were conducted on a square lattice 
with the coupling constant set to K = 0.1, otherwise all other parameters are the same as in d = 1. As in the previous 
case, the transition becomes sharper as one increases N. 



7 




FIG. 4: Mean field Z as a function of a for d = 2 for / x I lattices. I values are 7 (squares), 10 (pluses), 15 (diamonds), 20 
(circles). The coupling constant is K — 0.1, and results were averaged over more configurations, especially for low I. 



The phenomenon we present here is reminiscent of thermodynamic phase transitions in several features. First, 
the transition from unsynchronized state to a synchronized one breaks the original rotational symmetry of phases. 
Second, a phenomenon similar to the divergence of fluctuations close to the transition point can be observed also in 
this system, considering the standard deviation of the order parameter time-series (see Fig. pi). 
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FIG. 5: The standard deviation of the order parameter as a function of the a exponent for K =\. Symbols as in Fig. g. 

It is interesting to compare the behaviour of Z (Fig. ^) with the plot obtained looking at the fraction of oscillators 
locked at the mean frequency (Fig. used in the full locking approach. Although in the second case a discontinuity 
appears for finite N ^ one cannot distinguish between the transition region (with some unlocked oscillators but a strong 
mean field) and the region at high a (with developed incoherence and a low mean field). 
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FIG. 6: The fraction p of oscillators locked at the mean frequency for K = 1, N = 200. 

Another effect that was missed looking at full locking only is the dependence of the mean field on N . From Fig. |^ 
one can see that at low a Z increases with iV, since the oscillators rotate closer to each other (i.e, with smaller phase 
differences). We remark that this can happen even if the systems is not full locked. The effect on Z is the opposite 
for incoherence (high a). In this case, the decrease in the mean field is due to the fact that for high N , statistical 
fluctuations are reduced when averaging the uncorrelated phases over a larger number of the oscillators (Fig. ||). 
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FIG. 7: Enlargement of Fig. ^ Synchronization region. 




FIG. 8: Enlargement of Fig. g. Incoherence region. 
The dependence of Z over A'^ is also relevant when looking for a phase transition. In fact, it points to a discontinuity 
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in the plot Z = Z{a) for N going to infinity, as Figs. || and ^ suggest. In order to detect the phase transition, we 
perform simulations aimed at studying this effect. Results are plotted in Fig. ^. Considering as in Sec. II the 
convergence properties of the coupling term, and calling d the number of dimensions of the lattice, they can be 
interpreted as follows. 
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N 

FIG. 9: Size dependence of order parameter for a below (continuous line, a — 0.9) and above (dashed line, a = 1.1) the critical 
value a = d = 1. The coupling is K — 0.2 (triangles) and 0.1 (circles). 

For a > d the coupling term is bounded for N going to infinity. Hence, for K sufficiently small, the system is 
incoherent and the mean field approaches when N increases. For a < 1 the coupling term is unbounded for N going 
to infinity (i.e., may diverge in some regions of the phase space). Hence, for any (fixed) K, however small, the mean 
field asymptotically approaches 1 for N going to infinity. In mathematical terms: 

lim lim Z{a) = | i ^ ^ (8) 

This double limit gives a compact and meaningful result on the synchronization properties of a system in relation to 
the decay of the coupling signal. 

V. POPULATION OF OSCILLATORS WITH THERMAL NOISE 

The analogy with thermodynamic phase transitions can actually be developed further. In this Sec. we show that 
a similar transition takes place when the natural frequencies are equal, and randomness is introduced with thermal 
noise. The population of oscillators becomes then very close to a Heisenberg system with given temperature. 

We rewrite the natural frequencies in the form 

iU,{t) =LUo+^rit), (9) 

where ^i(t) is the noise term chosen from some distribution. We choose Cii^) a Gaussian distribution, such that 
£,i{t) — and £,i{t)^j{t') = 2DSij6{t — t'). It is clear that for every t the natural frequencies LOi have the same 
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distribution as in the above discussed model. However, the reahzation of this distribution changes at each instant of 
time. The equations of motion for the oscillators thus become: 



— sh 



^0 + ^ 2^ IF 



(10) 



Our simulations show that the phase transition from synchronized to unsynchronized state described in the previous 
Sec. takes place in this arrangement also (Figs. 10 and pl|) . Besides the transition, the size effects appear to remain 
valid in this case. For high _fC, synchronization does not depend on the relation between a and the lattice dimensions 
(Fig, p^ , triangles). However, for low coupling there is indication that the system synchronizes only if a < o? (Fig. 
O, circles). 
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FIG. 10: Average of the order parameter as a function of a for identical oscillators with noise in d = 1. Simulation results for 
7V=50, 100, 200, 500, and 1000 (symbols are squares, pluses, diamonds, circles, and down triangles respectively). The coupling 
is if = 1.0. 
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FIG. 11: Standard deviation of the order parameter as a function of a for identical oscillators with noise in d = 1. Symbols 
are the same as in Fig. |l^. 
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FIG. 12: Size dependence of order parameter for identical oscillators with noise for a below (continuous line) and above (dashed 
line) the critical value a = d = 1. The coupling is if = 1.0 (triangles) and 0.5 (circles). 
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VI. DISCUSSION 

We have considered synchronization of phase oscillators on a lattice, looking at critical levels of spatial decay in 
the interaction. Especially, since the coupling term considered is not normalized, we studied the effect of changing 
the population size on synchronization. We examined the system by investigating full locking and the mean field. 
The two approaches have appeared to be complementary. The criterion of full locking allowed to precisely define the 
boundary of complete synchronization for finite population sizes. It gives however a very strong condition, requiring 
all oscillators to be exactly locked. For small coupling, this condition is not useful. If the coupling constant is small, 
synchronization can still be enhanced, below a critical a value, enlarging the size of the population: but some of 
the oscillators remain unlocked. We thus studied the system from a difi^erent point of view, that is, looking at the 
behaviour of the mean field. In that case, a transition point is not strictly defined for finite population size. However, 
in the limit of infinite number of oscillators one can look, in analogy with thermodynamics, to the mean field as an 
order parameter and thus find critical values of the parameters where discontinuities appear. As one may expect, 
the value of the decay exponent equal to the number of lattice dimensions is then a good candidate for a transition 
point. At that value and below in fact, the coupling term is unbounded for an infinite size. That was supported by 
numerical simulations. As we pointed out, this gives a robust result for real systems: knowing only the number of 
lattice's dimensions and the decay in space of the coupling signal, one can predict if enlarging the size of the system 
eventually results in synchronization or not, even for arbitarily weak coupling constant. 

We finally considered a system of oscillators in which the diversity is given not by fixed natural frequencies, but by 
noise. The notion of full locking is not useful for this system, but the mean field approach can be carried out, and 
suggests the same features and critical point, at a decay exponent equal to the lattice's number of dimensions. Beside 
showing the robustness of the result, this last result is promising for an analytical treatment of the problem. In fact, 
the system with thermal noise is close to a Heisenberg spin-system. We remark that an analytical investigation would 
be important. In fact, simulations are very time consuming, and allowed us to have indications of the phenomenon in 
a relatively small parameter region. Especially, the analysis of higher dimensional lattices, as well as lower coupling 
strength, would be interesting. 
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